	implicit doubleprecision (a-h,o-z)
	dimension bf0(120,4),bf(120,4),bl(120,4),blf(120,4),bc(120,4)
	dimension r0(94),r(130)
	dimension bi0(94+1,25+1+1,4,7-1),bnl0(94+1,9+1,3,7),bm0(94+1,7,4)
	dimension curi0(4),curnl0(3),curm0(7)
	dimension ei0(94,25+1,4,7-1),enl0(94,9,3,7)
	dimension bi(130,25+1,4,7-1),bnl(130,9+1,3,7),bm(130,7,4)
	dimension ei(130,25,4,7-1),enl(130,9,3,7),em(130,7,4)
	dimension eibr(130,25,4,4),enlbr(130,9,3,4)
	dimension curi(25,4),curnl(9,4),curm(4)
	dimension eicr(130,25,4),enlcr(130,9,4),emcr(130,4)
	dimension ef0(130,35,4),ef(130*4,25+9*4+1)
	dimension h_matrix(62,62),b_victor(62)
	dimension g_victor(62,0:61),p_victor(62,0:61)
	dimension gthp(0:61),pthp(0:61),gtp(0:61)
	dimension alpha(0:61),r_lnbda(0:61)
	dimension curt0(35,0:1,4),curt(25+9*4+1,0:62),curl(25+9*4+1)
	dimension linkcb(35,35,4),noused(35,4)
	dimension qz(130,4)
	dimension inlm_over(35,4),lmover(25+9*4+1)
	integer is(4)
	character el*2,n_io*2,ppp*5
	character maglevel(7)*2,coe*3,signs*1
	data nr0,ns,nlv,ni,nnl,nn,nl / 94,4,7,25,9,5,4 /
	data (maglevel(l),l=1,7) / '06','10','12','14','15','16','17' /
	data curm0 / 1158.85d0, 1932.71d0, 2328.38d0, 2749.06d0,
	&                        2990.14d0, 3291.84d0, 3738.09d0 /
	data ( curi0(i),i=1,4), ncoei
	&   / -240.0d0,-120.0d0, 120.0d0, 240.0d0, 4 /
	data (curnl0(i),i=1,3),ncoenl
	&   /  80.0d0, 160.0d0, 240.0d0, 3 /
	data curip_limit,curnlp_limit,curmp_limit,
	&     curin_limit,curnln_limit,curmn_limit 
	&   / 240.0d0, 240.0d0, 4000.0d0,
	&    -240.0d0,-240.0d0, 1000.0d0 /
	common /cbf0/bf0
	common /cr0/r0
	common /cbi0/bi0
	common /cei0/ei0
	common /cbnl0/bnl0
	common /cenl0/enl0
	common /cbm0/bm0
	common /cbi/bi/cei/ei						   
	common /cbnl/bnl/cenl/enl
	common /ccuri0/curi0/ceibr/eibr
	common /ccurnl0/curnl0/cenlbr/enlbr
	common /ccurm0/curm0/cem/em
	common /cr/r/cbf/bf/cbc/bc

	open(unit=11,
	&file='n_it.dat',status='old')
	read(11,fmt='(2x,i2)') n_it
	close(11)
	n_io='00'
	if (n_it.eq.0) n_io='00'
	if (n_it.eq.1) n_io='01'
	if (n_it.eq.2) n_io='02'
	if (n_it.eq.3) n_io='03'
	if (n_it.eq.4) n_io='04'
	if (n_it.eq.5) n_io='05'
	if (n_it.eq.6) n_io='06'
	if (n_it.eq.7) n_io='07'
	if (n_it.eq.8) n_io='08'
	if (n_it.eq.9) n_io='09'
	if (n_it.eq.10) n_io='10'
	if (n_it.eq.11) n_io='11'
	if (n_it.eq.12) n_io='12'
	if (n_it.eq.13) n_io='13'
	if (n_it.eq.14) n_io='14'
	if (n_it.eq.15) n_io='15'
	open(unit=11,
	&file='saiifrev'//n_io//'.dat',status='old')
	read(11,fmt=
	&'(2x,a2,2x,a3,f6.1,2x,a3,f5.1,2x,a3,f7.2,a5,/,
	&2x,a4,f10.6,a3,2x,a3,f4.1,/)') 
	&el,ppp,a,ppp,z,ppp,w,ppp,
	&ppp,frad,ppp,ppp,h
	do n=1,ns
	read(11,fmt='(1x,i1)') is(n)
	read(11,fmt='(5(f10.3))') (bf0(i,n),i=1,120)
	end do
	close(11)
	do l=1,nlv
	open(unit=12,file='magmap\b'
	&//maglevel(l)//'a360p.dat',status='old')
	read(12,*)
	do n=1,ns
10	read(12,fmt='(5(f10.3))') (bm0(i,l,n),i=1,nr0+1)
	if (bm0(nr0+1,l,n).ne.90.0d0*(n-1)+45.0d0) goto 10
	end do
	close(12)
	end do
	do l=1,nlv-1
	do k=1,ncoei
	if (k.eq.1.or.k.eq.2) signs='n'
	if (k.eq.3.or.k.eq.4) signs='p'
	if (k.eq.1.or.k.eq.4) coe='240'
	if (k.eq.2.or.k.eq.3) coe='120'
	open (unit=13,file='mageff\bi'
	&//maglevel(l)//signs//coe//'.dat',status='old')
	do j=1,ni+1
	read (13,fmt='(1x,5(f10.3))') (bi0(i,j,k,l),i=1,nr0+1)
	end do
	close(13)
	end do
	end do

	open (unit=13,file='mageff\bi'
	&//maglevel(1)//'p'//'120'//'.dat',status='old')
	do j=1,ni+1+1
	read (13,fmt='(1x,5(f10.3))') (bi0(i,j,3,1),i=1,nr0+1)
	end do
	close(13)
	open (unit=13,file='mageff\bi'
	&//maglevel(2)//'n'//'240'//'.dat',status='old')
	do j=1,ni+1+1
	read (13,fmt='(1x,5(f10.3))') (bi0(i,j,1,2),i=1,nr0+1)
	end do
	close(13)
	open (unit=13,file='mageff\bi'
	&//maglevel(3)//'n'//'240'//'.dat',status='old')
	do j=1,ni+1+1
	read (13,fmt='(1x,5(f10.3))') (bi0(i,j,1,3),i=1,nr0+1)
	end do
	close(13)
	open (unit=13,file='mageff\bi'
	&//maglevel(5)//'n'//'240'//'.dat',status='old')
	do j=1,ni+1+1
	read (13,fmt='(1x,5(f10.3))') (bi0(i,j,1,5),i=1,nr0+1)
	end do
	close(13)

	signs='p'
	do l=1,nlv
	do k=1,ncoenl
	if (k.eq.1) coe='80'
	if (k.eq.2) coe='16'
	if (k.eq.3) coe='24'
	open (unit=13,file='mageff\bnl'
	&//maglevel(l)//signs//coe(1:2)//'.dat',status='old')
	do j=1,nn+nl+1
	read (13,fmt='(1x,5(f10.3))') (bnl0(i,j,k,l),i=1,nr0+1)
	end do
	close(13)
	end do
	end do

	do i=1,35
	r0(i)=84.0d0+2.0d0*dble(i-1)
	end do
	do i=36,71
	r0(i)=152.0d0+4.0d0*dble(i-35)
	end do
	do i=72,94
	r0(i)=296.0d0+2.0d0*dble(i-71)
	end do
	do l=1,nlv-1
	do k=1,ncoei
	do j=1,ni+1
	do i=1,nr0
	ei0(i,j,k,l)=(bi0(i,j+1,k,l)-bi0(i,j,k,l))/curi0(k)
	end do
	end do
	end do
	end do
	do l=1,nlv
	do k=1,ncoenl
	do j=1,nnl+1
	do i=1,nr0
	enl0(i,j,k,l)=(bnl0(i,j+1,k,l)-bnl0(i,j,k,l))/curnl0(k)
	end do
	end do
	end do
	end do

	do i=1,nr0
	eee=ei0(i,23,1,1)
	ei0(i,23,1,1)=ei0(i,24,1,1)
	ei0(i,24,1,1)=eee
	end do
	do j=19,ni
	do i=1,nr0
	ei0(i,j,3,1)=ei0(i,j+1,3,1)
	end do
	end do
	do j=24,ni
	do i=1,nr0
	ei0(i,j,1,2)=ei0(i,j+1,1,2)
	end do
	end do
	do j=11,ni
	do i=1,nr0
	ei0(i,j,1,3)=ei0(i,j+1,1,3)
	end do
	end do
	do j=23,ni
	do i=1,nr0
	ei0(i,j,1,5)=ei0(i,j+1,1,5)
	end do
	end do
	do i=1,nr0
	ei0(i,11,1,3)=ei0(i,11,1,3)-8.0d0/curi0(1)
	ei0(i,15,4,3)=ei0(i,15,3,3)
	ei0(i,18,4,3)=ei0(i,18,3,3)
	ei0(i, 4,3,4)=ei0(i, 4,3,4)/0.89044d0
	ei0(i, 4,4,4)=ei0(i, 4,4,4)/0.89044d0
	ei0(i,22,4,4)=ei0(i,22,4,4)*2.0d0
	end do
	do i=1,nr0
	enl0(i,9,3,3)=(enl0(i,9,1,3)*curnl0(1)+enl0(i,9,2,3)*curnl0(2))/
	&curnl0(3)
	end do

	nr=118
	do i=1,nr
	r(i)=104.0d0+2.0d0*dble(i-1)
	end do

	do n=1,ns
	do i=1,nr
	call larg1dr4bf(n,r(i),b)
	bf(i,n)=b
	bc(i,n)=0.0d0
	end do
	end do

	do l=1,nlv-1
	do k=1,ncoei
	do j=1,ni+1
	do i=1,nr
	call larg1dr4bi(l,k,j,r(i),b)
	bi(i,j,k,l)=b
	end do
	end do
	do j=1,ni
	do i=1,nr
	call larg1dr4ei(l,k,j,r(i),efi)
	ei(i,j,k,l)=efi
	end do
	end do
	end do
	end do
	do l=1,nlv
	do k=1,ncoenl
	do j=1,nnl+1
	do i=1,nr
	call larg1dr4bnl(l,k,j,r(i),b)
	bnl(i,j,k,l)=b
	end do
	end do
	do j=1,nnl
	do i=1,nr
	call larg1dr4enl(l,k,j,r(i),efnl)
	enl(i,j,k,l)=efnl
	end do
	end do
	end do
	end do
	do l=1,nlv
	do n=1,ns
	do i=1,nr
	call larg1dr4bm(n,l,r(i),b)
	bm(i,l,n)=b
	em(i,l,n)=bm(i,l,n)/curm0(l)
	end do
	end do
	end do

	do n=1,ns
	do k=1,ncoei
	do j=1,ni
	do i=1,nr
	call larg1dr4eibr(k,j,i,bf(i,n),efibr)
	eibr(i,j,k,n)=efibr
	end do
	end do
	end do
	do k=1,ncoenl
	do j=1,nnl
	do i=1,nr
	call larg1dr4enlbr(k,j,i,bf(i,n),efnlbr)
	enlbr(i,j,k,n)=efnlbr
	end do
	end do
	end do
	end do
	do n=1,ns
	do i=1,ni+nnl+1
	do j=1,ni+nnl+1
	linkcb(i,j,n)=0
	end do
	noused(i,n)=0
	end do
	end do
	linkcb(27,28,1)=1
	linkcb(29,30,1)=1
	linkcb(29,30,2)=1
	linkcb(29,30,4)=1
	do n=1,ns
	noused( 1,n)=1
	noused(25,n)=1
	end do

	do n=1,ns
	do i=1,nr
	qz(i,n)=1.0d0
	end do
	end do

	do n=1,ns
	do j=1,ni
	curi(j,n)=50.0d0
	end do
	do j=1,nnl
	curnl(j,n)=50.0d0
	end do
	curm(n)=2000.0d0
	end do

	do n=1,ns
	do j=1,ni
	if (noused(j,n).eq.1) then
	curi(j,n)=0.0d0
	end if
	end do
	end do
	do n=1,ns
	do i=ni+1,ni+nnl+1-2
	do j=i+1,ni+nnl+1-1
	if (linkcb(i,j,n).eq.1) then
	curnl(j,n)=curnl(i,n)
	end if
	end do
	end do
	end do

	accurate=1.0d-3
	relax_inc=2.0d1

	n_iterative=1

20	do n=1,ns
	do j=1,ni
	do i=1,nr
	call larg1dr4eicr(i,j,n,curi(j,n),eficr)
	eicr(i,j,n)=eficr
	end do
	end do
	do j=1,nnl
	do i=1,nr
	call larg1dr3enlcr(i,j,n,dabs(curnl(j,n)),efnlcr)
	enlcr(i,j,n)=efnlcr
	end do
	end do
	do i=1,nr
	call larg1dr4emcr(i,n,curm(n),efmcr)
	emcr(i,n)=efmcr
	end do
	end do

	do n=1,ns
	do j=1,ni
	do i=1,nr
	ef0(i,j,n)=eicr(i,j,n)
	end do
	curt0(j,0,n)=curi(j,n)
	end do
	do j=1,nnl
	do i=1,nr
	ef0(i,j+ni,n)=enlcr(i,j,n)
	end do
	curt0(j+ni,0,n)=curnl(j,n)
	end do
	do i=1,nr
	ef0(i,ni+nnl+1,n)=emcr(i,n)
	end do
	curt0(ni+nnl+1,0,n)=curm(n)
	end do

	do j=1,1+ni+nnl*ns
	do i=1,nr*ns
	ef(i,j)=0.0d0
	end do
	end do

	do n=1,ns
	do j=1,ni
	do i=1,nr
	ef(i+nr*(n-1),j)=ef0(i,j,n)
	end do
	curt(j,0)=curt0(j,0,n)
	end do
	do j=1+ni,nnl+ni
	do i=1,nr
	ef(i+nr*(n-1),j+nnl*(n-1))=ef0(i,j,n)
	end do
	curt(j+nnl*(n-1),0)=curt0(j,0,n)
	end do
	do i=1,nr
	ef(i+nr*(n-1),1+ni+nnl*ns)=ef0(i,1+ni+nnl,n)
	end do
	curt(1+ni+nnl*ns,0)=curt0(1+ni+nnl,0,n)
	end do

	lcb=0
	do n=ns,1,-1
	do ii=ni+nnl+1-2,ni+1,-1
	do jj=ni+nnl+1-1,ii+1,-1
	if (linkcb(ii,jj,n).eq.1) then
	do i=1+nr*(n-1),nr+nr*(n-1)
	ef(i,ii+nnl*(n-1))=ef(i,ii+nnl*(n-1))+ef(i,jj+nnl*(n-1))
	end do
	do j=jj+nnl*(n-1),62-1-lcb
	do i=1,nr*ns
	ef(i,j)=ef(i,j+1)
	end do
	curt(j,0)=curt(j+1,0)
	end do
	lcb=lcb+1
	end if
	end do
	end do
	end do

	nud=0
	n=1
	do jj=ni,1,-1
	if (noused(jj,n).eq.1) then
	do j=jj,62-1-lcb-nud
	do i=1,nr*ns
	ef(i,j)=ef(i,j+1)
	end do
	curt(j,0)=curt(j+1,0)
	end do
	nud=nud+1
	end if
	end do

	lcb0=lcb
	nud0=nud
	
	do i=1,62
	do j=1,62
	h_matrix(i,j)=0.0d0
	end do
	b_victor(i)=0.0d0
	end do
	do i=1,62-lcb0-nud0
	do j=1,62-lcb0-nud0
	do ii=1,nr*ns
	h_matrix(i,j)=h_matrix(i,j)+2.0d0*ef(ii,i)*ef(ii,j)
	&*qz(mod(int4(ii-1),int4(nr))+1,(ii-1)/nr+1)
	end do
	end do
	do ii=1,nr*ns
	b_victor(i)=b_victor(i)+2.0d0*ef(ii,i)
	&*(0.0d0-bf(mod(int4(ii-1),int4(nr))+1,(ii-1)/nr+1))
	&*qz(mod(int4(ii-1),int4(nr))+1,(ii-1)/nr+1)
	end do
	end do

	do k=1-1,62-1
	do i=1,62
	g_victor(i,k)=0.0d0
	p_victor(i,k)=0.0d0
	end do
	end do

	do k=1-1,62-1-lcb0-nud0

	do kk=1-1,62-1
	gthp(kk)=0.0d0
	gtp(kk)=0.0d0
	pthp(kk)=0.0d0
	end do

	do i=1,62-lcb0-nud0
	do j=1,62-lcb0-nud0
	g_victor(i,k)=g_victor(i,k)+h_matrix(i,j)*curt(j,k)
	end do
	g_victor(i,k)=g_victor(i,k)+b_victor(i)
	end do
	if (k.eq.0) then
	do i=1,62-lcb0-nud0
	p_victor(i,k)=0.0d0-g_victor(i,k)
	end do
	end if
	if (k.gt.0) then
	do kk=0,k-1
	do i=1,62-lcb0-nud0
	do j=1,62-lcb0-nud0
	gthp(kk)=gthp(kk)+g_victor(i,k)*h_matrix(i,j)*p_victor(j,kk)
	pthp(kk)=pthp(kk)+p_victor(i,kk)*h_matrix(i,j)*p_victor(j,kk)
	end do
	end do
	alpha(kk)=gthp(kk)/pthp(kk)
	do i=1,62-lcb0-nud0
	p_victor(i,k)=p_victor(i,k)+alpha(kk)*p_victor(i,kk)
	end do
	end do
	do i=1,62-lcb0-nud0
	p_victor(i,k)=p_victor(i,k)-g_victor(i,k)
	end do
	end if
	do i=1,62-lcb0-nud0
	gtp(k)=gtp(k)+g_victor(i,k)*p_victor(i,k)
	do j=1,62-lcb0-nud0
	pthp(k)=pthp(k)+p_victor(i,k)*h_matrix(i,j)*p_victor(j,k)
	end do
	end do
	r_lnbda(k)=0.0d0-gtp(k)/pthp(k)
	do i=1,62-lcb0-nud0
	curt(i,k+1)=curt(i,k)+r_lnbda(k)*p_victor(i,k)
	end do

	end do
	
	n=1
	do jj=1,ni
	if (noused(jj,n).eq.1) then
	do j=62-lcb-nud,jj,-1
	curt(j+1,62-lcb0-nud0)=curt(j,62-lcb0-nud0)
	end do
	curt(jj,62-lcb0-nud0)=0.0d0
	do nn=1,ns
	curt0(jj,0,nn)=0.0d0
	end do
	nud=nud-1
	end if
	end do

	do n=1,ns
	do ii=ni+1,ni+nnl+1-2
	do jj=ii+1,ni+nnl+1-1
	if (linkcb(ii,jj,n).eq.1) then
	do j=62-lcb,jj+nnl*(n-1),-1
	curt(j+1,62-lcb0-nud0)=curt(j,62-lcb0-nud0)
	end do
	curt(jj+nnl*(n-1),62-lcb0-nud0)=curt(ii+nnl*(n-1),62-lcb0-nud0)
	curt0(jj,0,n)=curt(ii+nnl*(n-1),62-lcb0-nud0)
	lcb=lcb-1
	end if
	end do
	end do
	end do

	do n=1,ns
	do j=1,ni
	curt0(j,1,n)=curt(j,62-lcb0-nud0)
	end do
	do j=1,nnl*ns
	curt0(ni+mod(int4(j-1),int4(nnl))+1,1,(j-1)/nnl+1)=
	&curt(j+ni,62-lcb0-nud0)
	end do
	curt0(ni+nnl+1,1,n)=curt(ni+nnl*ns+1,62-lcb0-nud0)
	end do

	do n=1,ns
	do j=1,ni
	curi(j,n)=curi(j,n)+(curt0(j,1,n)-curi(j,n))/relax_inc
	end do
	do j=1,nnl
	curnl(j,n)=curnl(j,n)+(curt0(j+ni,1,n)-curnl(j,n))/relax_inc
	end do
	curm(n)=curm(n)+(curt0(ni+nnl+1,1,n)-curm(n))/relax_inc
	end do

	do n=1,ns
	do j=1,ni+nnl+1
	if ( dabs(curt0(j,1,n)-curt0(j,0,n)).gt.accurate ) then
	write (6,fmt='(3x,a9,i2,3x,a7,i3,3x,a10,f12.5,3x,a12,i7)') 
	&'n_sector=',n,'n_coil=',j,
	&'d_current=',(curt0(j,1,n)-curt0(j,0,n)),
	&'n_iterative=',n_iterative
	n_iterative=n_iterative+1
	goto 20
	end if
	end do
	end do

	do n=1,ns
	do i=1,nr
	do j=1,ni+nnl+1
	bc(i,n)=bc(i,n)+ef0(i,j,n)*curt0(j,1,n)
	end do
	end do
	end do

	open (unit=14,
	&file='result00.dat',status='unknown')
	do j=1,ni
	write(14,fmt='(4(3x,a5,i2,a1,i1,a2,f9.3))') 
	&('curi(',j,',',n,')=',curi(j,n),n=1,ns)
	end do
	do j=1,nnl
	write(14,fmt='(4(3x,a6,i1,a1,i1,a2,f9.3))')
	&('curnl(',j,',',n,')=',curnl(j,n),n=1,ns)
	end do
	write(14,fmt='(4(6x,a5,i1,a2,f9.3),/)')
	&('curm(',n,')=',curm(n),n=1,ns)
	do n=1,ns
	do i=1,nr
	write(14,fmt='(7x,a2,i3,a2,f9.3,3(3x,a4,i3,a1,i1,a2,f9.3))')
	&'r(',i,')=',r(i),
	&'bfc(',i,',',n,')=',bc(i,n),'bfi(',i,',',n,')=',bf(i,n),
	&'dtb(',i,',',n,')=',bc(i,n)-bf(i,n)
	end do
	end do
	close(14)

	call paint_bf()

	accurate=1.0d-3
	relax_inc=2.0d1

	do n=1,ns
	do j=1,1+ni+nnl
	inlm_over(j,n)=0
	end do
	end do

	n_overlimit=0

30	curp_max=0.0d0
	curn_max=0.0d0
	j_p=0
	n_p=0
	n=1
	do j=1,ni
	if (curi(j,n).gt.curip_limit.and.curi(j,n).gt.curp_max) then
	curp_max=curi(j,n)
	j_p=j
	n_p=n
	end if
	end do
	do n=1,ns
	do j=1,nnl
	if (curnl(j,n).gt.curnlp_limit.and.curnl(j,n).gt.curp_max) then
	curp_max=curnl(j,n)
	j_p=j+ni
	n_p=n
	end if
	end do
	end do
	j_n=0
	n_n=0
	n=1
	do j=1,ni
	if (curi(j,n).lt.curin_limit.and.curi(j,n).lt.curn_max) then
	curn_max=curi(j,n)
	j_n=j
	n_n=n
	end if
	end do
	do n=1,ns
	do j=1,nnl
	if (curnl(j,n).lt.curnln_limit.and.curnl(j,n).lt.curn_max) then
	curn_max=curnl(j,n)
	j_n=j+ni
	n_n=n
	end if
	end do
	end do

	if ( (curp_max-curip_limit).ge.(curin_limit-curn_max) ) then
	if ( j_p.ne.0.and.j_n.ne.0 ) then
	j_over=j_p
	n_over=n_p
	end if
	else
	if ( j_p.ne.0.and.j_n.ne.0 ) then
	j_over=j_n
	n_over=n_n
	end if
	end if
	if ( j_p.ne.0.and.j_n.eq.0 ) then
	j_over=j_p
	n_over=n_p
	end if
	if ( j_p.eq.0.and.j_n.ne.0 ) then
	j_over=j_n
	n_over=n_n
	end if
	n=1
	if ( curm(n).gt.curmp_limit.or.curm(n).lt.curmn_limit ) then
	if ( j_p.eq.0.and.j_n.eq.0 ) then
	j_over=1+ni+nnl
	n_over=n
	end if
	else
	if ( j_p.eq.0.and.j_n.eq.0 ) then
	j_over=0
	n_over=0
	end if
	end if

	if ( j_over.eq.0.and.n_over.eq.0 ) then
	goto 50
	else
	inlm_over(j_over,n_over)=1
	n_overlimit=n_overlimit+1
	end if

	n=1
	do j=1,ni
	if ( inlm_over(j,n).eq.1.and.curi(j,n).gt.curip_limit) then
	do nn=1,ns
	curi(j,nn)=curip_limit
	end do
	end if
	if ( inlm_over(j,n).eq.1.and.curi(j,n).lt.curin_limit) then
	do nn=1,ns
	curi(j,nn)=curin_limit
	end do
	end if
	end do
	do n=1,ns
	do j=1+ni,nnl+ni
	if ( inlm_over(j,n).eq.1.and.curnl(j-ni,n).gt.curnlp_limit) then
	curnl(j-ni,n)=curnlp_limit
	end if
	if ( inlm_over(j,n).eq.1.and.curnl(j-ni,n).lt.curnln_limit) then
	curnl(j-ni,n)=curnln_limit
	end if
	end do
	end do
	n=1
	j=1+ni+nnl
	if ( inlm_over(j,n).eq.1.and.curm(n).gt.curmp_limit) then
	do nn=1,ns
	curm(nn)=curmp_limit
	end do
	end if
	if ( inlm_over(j,n).eq.1.and.curm(n).lt.curmn_limit) then
	do nn=1,ns
	curm(nn)=curmn_limit
	end do
	end if

	do n=1,ns
	do j=1,ni
	if (noused(j,n).eq.1) then
	curi(j,n)=0.0d0
	end if
	end do
	end do
	do n=1,ns
	do i=ni+1,ni+nnl+1-2
	do j=i+1,ni+nnl+1-1
	if (linkcb(i,j,n).eq.1) then
	curnl(j,n)=curnl(i,n)
	end if
	end do
	end do
	end do

	n_iterative=1

40	do n=1,ns
	do j=1,ni
	do i=1,nr
	call larg1dr4eicr(i,j,n,curi(j,n),eficr)
	eicr(i,j,n)=eficr
	end do
	end do
	do j=1,nnl
	do i=1,nr
	call larg1dr3enlcr(i,j,n,dabs(curnl(j,n)),efnlcr)
	enlcr(i,j,n)=efnlcr
	end do
	end do
	do i=1,nr
	call larg1dr4emcr(i,n,curm(n),efmcr)
	emcr(i,n)=efmcr
	end do
	end do

	do n=1,ns
	do j=1,ni
	do i=1,nr
	ef0(i,j,n)=eicr(i,j,n)
	end do
	curt0(j,0,n)=curi(j,n)
	end do
	do j=1,nnl
	do i=1,nr
	ef0(i,j+ni,n)=enlcr(i,j,n)
	end do
	curt0(j+ni,0,n)=curnl(j,n)
	end do
	do i=1,nr
	ef0(i,ni+nnl+1,n)=emcr(i,n)
	end do
	curt0(ni+nnl+1,0,n)=curm(n)
	end do

	do j=1,1+ni+nnl*ns
	do i=1,nr*ns
	ef(i,j)=0.0d0
	end do
	end do

	do n=1,ns
	do j=1,ni
	do i=1,nr
	ef(i+nr*(n-1),j)=ef0(i,j,n)
	end do
	curt(j,0)=curt0(j,0,n)
	end do
	do j=1+ni,nnl+ni
	do i=1,nr
	ef(i+nr*(n-1),j+nnl*(n-1))=ef0(i,j,n)
	end do
	curt(j+nnl*(n-1),0)=curt0(j,0,n)
	end do
	do i=1,nr
	ef(i+nr*(n-1),1+ni+nnl*ns)=ef0(i,1+ni+nnl,n)
	end do
	curt(1+ni+nnl*ns,0)=curt0(1+ni+nnl,0,n)
	end do

	n=1
	do j=1,ni
	lmover(j)=inlm_over(j,n)
	end do
	do n=1,ns
	do j=1+ni,nnl+ni
	lmover(j+(n-1)*nnl)=inlm_over(j,n)
	end do
	end do
	n=1
	lmover(ni+ns*nnl+1)=inlm_over(ni+nnl+1,n)

	lcb=0
	do n=ns,1,-1
	do ii=ni+nnl+1-2,ni+1,-1
	do jj=ni+nnl+1-1,ii+1,-1
	if (linkcb(ii,jj,n).eq.1) then
	do i=1+nr*(n-1),nr+nr*(n-1)
	ef(i,ii+nnl*(n-1))=ef(i,ii+nnl*(n-1))+ef(i,jj+nnl*(n-1))
	end do
	do j=jj+nnl*(n-1),62-1-lcb
	do i=1,nr*ns
	ef(i,j)=ef(i,j+1)
	end do
	curt(j,0)=curt(j+1,0)
	lmover(j)=lmover(j+1)
	end do
	lcb=lcb+1
	end if
	end do
	end do
	end do

	nud=0
	n=1
	do jj=ni,1,-1
	if (noused(jj,n).eq.1) then
	do j=jj,62-1-lcb-nud
	do i=1,nr*ns
	ef(i,j)=ef(i,j+1)
	end do
	curt(j,0)=curt(j+1,0)
	lmover(j)=lmover(j+1)
	end do
	nud=nud+1
	end if
	end do

	do n=1,ns
	do i=1,nr
	bl(i,n)=0.0d0
	end do
	end do

	lor=0
	jj=ni+ns*nnl+1-lcb-nud
	if (lmover(jj).eq.1) then
	do i=1,ns*nr
	bl(mod(int4(i-1),int4(nr))+1,(i-1)/nr+1)=
	&bl(mod(int4(i-1),int4(nr))+1,(i-1)/nr+1)+ef(i,jj)*curt(jj,0)
	end do
	curl(jj)=curt(jj,0)
	lor=lor+1
	end if
	do jj=ni+ns*nnl+1-1-lcb-nud,1,-1
	if (lmover(jj).eq.1) then
	do i=1,ns*nr
	bl(mod(int4(i-1),int4(nr))+1,(i-1)/nr+1)=
	&bl(mod(int4(i-1),int4(nr))+1,(i-1)/nr+1)+ef(i,jj)*curt(jj,0)
	end do
	curl(jj)=curt(jj,0)
	do j=jj,62-1-lcb-nud-lor
	do i=1,nr*ns
	ef(i,j)=ef(i,j+1)
	end do
	curt(j,0)=curt(j+1,0)
	end do
	lor=lor+1
	end if
	end do

	do n=1,ns
	do i=1,nr
	blf(i,n)=bf(i,n)-bl(i,n)
	end do
	end do

	lcb0=lcb
	nud0=nud
	lor0=lor
	
	do i=1,62
	do j=1,62
	h_matrix(i,j)=0.0d0
	end do
	b_victor(i)=0.0d0
	end do
	do i=1,62-lcb0-nud0-lor0
	do j=1,62-lcb0-nud0-lor0
	do ii=1,nr*ns
	h_matrix(i,j)=h_matrix(i,j)+2.0d0*ef(ii,i)*ef(ii,j)
	&*qz(mod(int4(ii-1),int4(nr))+1,(ii-1)/nr+1)
	end do
	end do
	do ii=1,nr*ns
	b_victor(i)=b_victor(i)+2.0d0*ef(ii,i)
	&*(0.0d0-blf(mod(int4(ii-1),int4(nr))+1,(ii-1)/nr+1))
	&*qz(mod(int4(ii-1),int4(nr))+1,(ii-1)/nr+1)
	end do
	end do

	do k=1-1,62-1
	do i=1,62
	g_victor(i,k)=0.0d0
	p_victor(i,k)=0.0d0
	end do
	end do

	do k=1-1,62-1-lcb0-nud0-lor0

	do kk=1-1,62-1
	gthp(kk)=0.0d0
	gtp(kk)=0.0d0
	pthp(kk)=0.0d0
	end do

	do i=1,62-lcb0-nud0-lor0
	do j=1,62-lcb0-nud0-lor0
	g_victor(i,k)=g_victor(i,k)+h_matrix(i,j)*curt(j,k)
	end do
	g_victor(i,k)=g_victor(i,k)+b_victor(i)
	end do
	if (k.eq.0) then
	do i=1,62-lcb0-nud0-lor0
	p_victor(i,k)=0.0d0-g_victor(i,k)
	end do
	end if
	if (k.gt.0) then
	do kk=0,k-1
	do i=1,62-lcb0-nud0-lor0
	do j=1,62-lcb0-nud0-lor0
	gthp(kk)=gthp(kk)+g_victor(i,k)*h_matrix(i,j)*p_victor(j,kk)
	pthp(kk)=pthp(kk)+p_victor(i,kk)*h_matrix(i,j)*p_victor(j,kk)
	end do
	end do
	alpha(kk)=gthp(kk)/pthp(kk)
	do i=1,62-lcb0-nud0-lor0
	p_victor(i,k)=p_victor(i,k)+alpha(kk)*p_victor(i,kk)
	end do
	end do
	do i=1,62-lcb0-nud0-lor0
	p_victor(i,k)=p_victor(i,k)-g_victor(i,k)
	end do
	end if
	do i=1,62-lcb0-nud0-lor0
	gtp(k)=gtp(k)+g_victor(i,k)*p_victor(i,k)
	do j=1,62-lcb0-nud0-lor0
	pthp(k)=pthp(k)+p_victor(i,k)*h_matrix(i,j)*p_victor(j,k)
	end do
	end do
	r_lnbda(k)=0.0d0-gtp(k)/pthp(k)
	do i=1,62-lcb0-nud0-lor0
	curt(i,k+1)=curt(i,k)+r_lnbda(k)*p_victor(i,k)
	end do

	end do

	do jj=1,ni+ns*nnl+1-1-lcb-nud
	if (lmover(jj).eq.1) then
	do j=62-lcb-nud-lor,jj,-1
	curt(j+1,62-lcb0-nud0-lor0)=curt(j,62-lcb0-nud0-lor0)
	end do
	curt(jj,62-lcb0-nud0-lor0)=curl(jj)
	lor=lor-1
	end if
	end do
	jj=ni+ns*nnl+1-lcb-nud
	if (lmover(jj).eq.1) then
	curt(jj,62-lcb0-nud0-lor0)=curl(jj)
	lor=lor-1
	end if
	
	n=1
	do jj=1,ni
	if (noused(jj,n).eq.1) then
	do j=62-lcb-nud,jj,-1
	curt(j+1,62-lcb0-nud0-lor0)=curt(j,62-lcb0-nud0-lor0)
	end do
	curt(jj,62-lcb0-nud0-lor0)=0.0d0
	do nn=1,ns
	curt0(jj,0,nn)=0.0d0
	end do
	nud=nud-1
	end if
	end do

	do n=1,ns
	do ii=ni+1,ni+nnl+1-2
	do jj=ii+1,ni+nnl+1-1
	if (linkcb(ii,jj,n).eq.1) then
	do j=62-lcb,jj+nnl*(n-1),-1
	curt(j+1,62-lcb0-nud0-lor0)=curt(j,62-lcb0-nud0-lor0)
	end do
	curt(jj+nnl*(n-1),62-lcb0-nud0-lor0)=
	&curt(ii+nnl*(n-1),62-lcb0-nud0-lor0)
	curt0(jj,0,n)=curt(ii+nnl*(n-1),62-lcb0-nud0-lor0)
	lcb=lcb-1
	end if
	end do
	end do
	end do

	do n=1,ns
	do j=1,ni
	curt0(j,1,n)=curt(j,62-lcb0-nud0-lor0)
	end do
	do j=1,nnl*ns
	curt0(ni+mod(int4(j-1),int4(nnl))+1,1,(j-1)/nnl+1)=
	&curt(j+ni,62-lcb0-nud0-lor0)
	end do
	curt0(ni+nnl+1,1,n)=curt(ni+nnl*ns+1,62-lcb0-nud0-lor0)
	end do

	do n=1,ns
	do j=1,ni
c	curi(j,n)=curt0(j,ni+nnl+1,n)
	curi(j,n)=curi(j,n)+(curt0(j,1,n)-curi(j,n))/relax_inc
	end do
	do j=1,nnl
c	curnl(j,n)=curt0(j+ni,ni+nnl+1,n)
	curnl(j,n)=curnl(j,n)+(curt0(j+ni,1,n)-curnl(j,n))/relax_inc
	end do
c	curm(n)=curt0(ni+nnl+1,ni+nnl+1,n)
	curm(n)=curm(n)+(curt0(ni+nnl+1,1,n)-curm(n))/relax_inc
	end do

	do n=1,ns
	do j=1,ni+nnl+1
	if ( dabs(curt0(j,1,n)-curt0(j,0,n)).gt.accurate ) then
	write (6,
	&fmt='(3x,a9,i2,3x,a7,i3,3x,a10,f12.5,3x,a12,i7,3x,a12,i3)') 
	&'n_sector=',n,'n_coil=',j,
	&'d_current=',(curt0(j,1,n)-curt0(j,0,n)),
	&'n_iterative=',n_iterative,'n_overlimit=',n_overlimit
	n_iterative=n_iterative+1
	goto 40
	end if
	end do
	end do

	goto 30

50	do n=1,ns
	do i=1,nr
	bc(i,n)=0.0d0
	end do
	end do

	do n=1,ns
	do i=1,nr
	do j=1,ni+nnl+1
	bc(i,n)=bc(i,n)+ef0(i,j,n)*curt0(j,1,n)
	end do
	end do
	end do
	
	open (unit=15,
	&file='result'//n_io//'.dat',status='unknown')
	do j=1,ni
	write(15,fmt='(4(3x,a5,i2,a1,i1,a2,f9.3))') 
	&('curi(',j,',',n,')=',curi(j,n),n=1,ns)
	end do
	do j=1,nnl
	write(15,fmt='(4(3x,a6,i1,a1,i1,a2,f9.3))')
	&('curnl(',j,',',n,')=',curnl(j,n),n=1,ns)
	end do
	write(15,fmt='(4(6x,a5,i1,a2,f9.3),/)')
	&('curm(',n,')=',curm(n),n=1,ns)
	do n=1,ns
	do i=1,nr
	write(15,fmt='(7x,a2,i3,a2,f9.3,3(3x,a4,i3,a1,i1,a2,f9.3))')
	&'r(',i,')=',r(i),
	&'bfc(',i,',',n,')=',bc(i,n),'bfi(',i,',',n,')=',bf(i,n),
	&'dtb(',i,',',n,')=',bc(i,n)-bf(i,n)
	end do
	end do
	close(15)

	call paint_bf()

	stop
	end

	subroutine larg1dr4bf(n,x,b)
	implicit real(8) (a-h,o-z)
	dimension rb0(120),bf0(120,4)
	common /cbf0/bf0
	dimension x0(120),xi(4)
	do i=1,120
	rb0(i)=100.0d0+2.0d0*dble(i-1)
	end do
	ix_max=120
	do i=1,ix_max
	x0(i)=rb0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	b=0.0d0
	do i=1,4
	b=b+xi(i)*bf0(ix+i,n)
	end do
	return
	end

	subroutine larg1dr4bi(ll,kk,jj,x,b)
	implicit real(8) (a-h,o-z)
	dimension r0(94),bi0(94+1,25+1+1,4,7-1)
	common /cr0/r0/cbi0/bi0
	dimension x0(94),xi(4)
	ix_max=94
	do i=1,ix_max
	x0(i)=r0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	b=0.0d0
	do i=1,4
	b=b+xi(i)*bi0(ix+i,jj,kk,ll)
	end do
	return
	end

	subroutine larg1dr4ei(ll,kk,jj,x,e)
	implicit real(8) (a-h,o-z)
	dimension r0(94),ei0(94,25+1,4,7-1)
	common /cr0/r0/cei0/ei0
	dimension x0(94),xi(4)
	ix_max=94
	do i=1,ix_max
	x0(i)=r0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*ei0(ix+i,jj,kk,ll)
	end do
	return
	end

	subroutine larg1dr4bnl(ll,kk,jj,x,b)
	implicit real(8) (a-h,o-z)
	dimension r0(94),bnl0(94+1,9+1,3,7)
	common /cr0/r0/cbnl0/bnl0
	dimension x0(94),xi(4)
	ix_max=94
	do i=1,ix_max
	x0(i)=r0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	b=0.0d0
	do i=1,4
	b=b+xi(i)*bnl0(ix+i,jj,kk,ll)
	end do
	return
	end

	subroutine larg1dr4enl(ll,kk,jj,x,e)
	implicit real(8) (a-h,o-z)
	dimension r0(94),enl0(94,9,3,7)
	common /cr0/r0/cenl0/enl0
	dimension x0(94),xi(4)
	ix_max=94
	do i=1,ix_max
	x0(i)=r0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*enl0(ix+i,jj,kk,ll)
	end do
	return
	end

	subroutine larg1dr4bm(nn,ll,x,b)
	implicit real(8) (a-h,o-z)
	dimension r0(94),bm0(94+1,7,4)
	common /cr0/r0/cbm0/bm0
	dimension x0(94),xi(4)
	ix_max=94
	do i=1,ix_max
	x0(i)=r0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	b=0.0d0
	do i=1,4
	b=b+xi(i)*bm0(ix+i,ll,nn)
	end do
	return
	end

	subroutine larg1dr4eibr(kk,jj,ii,x,e)
	implicit real(8) (a-h,o-z)
	dimension bi(130,25+1,4,7-1),ei(130,25,4,7-1)
	common /cbi/bi/cei/ei						   
	dimension x0(7-1),xi(4)
	ix_max=7-1
	do i=1,ix_max
	x0(i)=bi(ii,jj+1,kk,i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*ei(ii,jj,kk,ix+i)
	end do
	return
	end

	subroutine larg1dr4enlbr(kk,jj,ii,x,e)
	implicit real(8) (a-h,o-z)
	dimension bnl(130,9+1,3,7),enl(130,9,3,7)
	common /cbnl/bnl/cenl/enl
	dimension x0(7),xi(4)
	ix_max=7
	do i=1,ix_max
	x0(i)=bnl(ii,jj+1,kk,i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*enl(ii,jj,kk,ix+i)
	end do
	return
	end

	subroutine larg1dr4embr(nn,ii,x,e)
	implicit real(8) (a-h,o-z)
	dimension bm(130,7,4),em(130,7,4)
	common /cbm/bm/cem/em
	dimension x0(7),xi(4)
	ix_max=7
	do i=1,ix_max
	x0(i)=bm(ii,i,nn)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*em(ii,ix+i,nn)
	end do
	return
	end

	subroutine larg1dr4eicr(ii,jj,nn,x,e)
	implicit real(8) (a-h,o-z)
	dimension curi0(4),eibr(130,25,4,4)
	common /ccuri0/curi0/ceibr/eibr
	dimension x0(4),xi(4)
	ix_max=4
	do i=1,ix_max
	x0(i)=curi0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*eibr(ii,jj,ix+i,nn)
	end do
	return
	end

	subroutine larg1dr3enlcr(ii,jj,nn,x,e)
	implicit real(8) (a-h,o-z)
	dimension curnl0(3),enlbr(130,9,3,4)
	common /ccurnl0/curnl0/cenlbr/enlbr
	dimension x0(3),xi(3)
	ix_max=3
	do i=1,ix_max
	x0(i)=curnl0(i)
	end do
	ix=0
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	xi(1)=(x-x2)*(x-x3)/((x1-x2)*(x1-x3))
	xi(2)=(x-x1)*(x-x3)/((x2-x1)*(x2-x3))
	xi(3)=(x-x1)*(x-x2)/((x3-x1)*(x3-x2))
	e=0.0d0
	do i=1,3
	e=e+xi(i)*enlbr(ii,jj,ix+i,nn)
	end do
	return
	end

	subroutine larg1dr4emcr(ii,nn,x,e)
	implicit real(8) (a-h,o-z)
	dimension curm0(7),em(130,7,4)
	common /ccurm0/curm0/cem/em
	dimension x0(7),xi(4)
	ix_max=7
	do i=1,ix_max
	x0(i)=curm0(i)
	end do
	do i=1,ix_max-1
	if (x0(i).le.x.and.x.lt.x0(i+1)) ix=i-2
	end do
	if (x.le.x0(1+1)) ix=0
	if (x.ge.x0(ix_max-1)) ix=ix_max-4
	x1=x0(ix+1)
	x2=x0(ix+2)
	x3=x0(ix+3)
	x4=x0(ix+4)
	xi(1)=(x-x2)*(x-x3)*(x-x4)/((x1-x2)*(x1-x3)*(x1-x4))
	xi(2)=(x-x1)*(x-x3)*(x-x4)/((x2-x1)*(x2-x3)*(x2-x4))
	xi(3)=(x-x1)*(x-x2)*(x-x4)/((x3-x1)*(x3-x2)*(x3-x4))
	xi(4)=(x-x1)*(x-x2)*(x-x3)/((x4-x1)*(x4-x2)*(x4-x3))
	e=0.0d0
	do i=1,4
	e=e+xi(i)*em(ii,ix+i,nn)
	end do
	return
	end

	subroutine lag4(n,xi,yi,x,y,io)
	implicit real(8) (a-h,o-z)
	dimension xi(n),yi(n),p(4)
	integer n,io
	do i=1,n
	if (x.eq.xi(i)) then
	y=yi(i)
	return
	end if
	end do
	io=0
	k =0
	if (x.lt.xi(1)) io=-1
	if (x.gt.xi(n)) io= 1
	do i=1,n-1
	if (x.gt.xi(i).and.x.lt.xi(i+1)) k=i
	end do
	kk=k-1
	if (kk.lt.0) kk=1
	if (kk.gt.n-3) kk=n-3
	if (io.eq.1) kk=n-3
	sum=0.0d0
	do i=1,4
	p(i)=1.0d0
	do j=1,4
	if (i.eq.j) goto 10
	p(i)=p(i)*(x-xi(j+kk-1))/(xi(i+kk-1)-xi(j+kk-1))
10	end do
	sum=sum+p(i)*yi(i+kk-1)
	end do
	y=sum
	return
	end

	subroutine paint_ef()
	use msflib
	implicit doubleprecision (a-h,o-z)
	dimension r0(94),r(130),ei0(94,25+1,4,7-1),enl0(94,9,3,7)
c	dimension bi(130,25+1,4,7-1),bnl(130,9+1,3,7),bm(130,7,4)
	dimension ei(130,25,4,7-1),enl(130,9,3,7)!,em(130,7,4)
	common /cr0/r0/cei0/ei0/cenl0/enl0
	common /cr/r/cei/ei/cenl/enl

	LOGICAL               statusmode
	TYPE (windowconfig)   myscreen

	INTEGER(2)            status
	INTEGER(2)            xwidth, yheight, cols, rows
      
	INTEGER(4)            rgbcolor, oldcolor, rgbbkcolor, oldbkcolor
	TYPE (wxycoord)       wxy
	TYPE (rccoord)        curpos

	integer(4)            result

!     Set the screen to the best resolution and maximum number of 
!     available colors.
	myscreen.numxpixels   = -1
	myscreen.numypixels   = -1
	myscreen.numtextcols  = -1
	myscreen.numtextrows  = -1
	myscreen.numcolors    = -1
	myscreen.fontsize = -1
	myscreen.title = " "C
	statusmode = SETWINDOWCONFIG(myscreen)
	IF(.NOT. statusmode) statusmode = SETWINDOWCONFIG(myscreen)

	statusmode = GETWINDOWCONFIG( myscreen )
	xwidth  = myscreen.numxpixels
	yheight = myscreen.numypixels
	cols    = myscreen.numtextcols
	rows    = myscreen.numtextrows

c	rgbbkcolor = #ffffff
	rgbbkcolor = #000000
	oldbkcolor = SETBKCOLORRGB( rgbbkcolor )
	  
	CALL CLEARSCREEN( $GCLEARSCREEN )

c	CALL SETVIEWPORT( INT2(0), INT2(0), xwidth-1, yheight-1 )
	CALL SETVIEWPORT( INT2(30), INT2(30), xwidth-1-30, yheight-1-30 )
	CALL SETTEXTWINDOW( INT2(1), INT2(1), rows, cols )
	x_min=r0(1)
	x_max=r0(94)
	y_min=-1.0d0
	y_max=1.0d0
	status = SETWINDOW( .true., x_min, y_min, x_max, y_max)

c	rgbcolor = #000000
	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )  
	status = RECTANGLE( $GBORDER,INT2(0),INT2(0),xwidth-1,yheight-1 )
!
!     Plot the points.
!
	
	do l=1,7-1
	do k=1,4
	do j=1,25
	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )
	status = RECTANGLE_W ( $GBORDER, x_min, y_min, x_max, y_max )
	CALL MOVETO_W ( x_min, (y_min+y_max)/2.0d0, wxy )
	result = LINETO_W ( x_max, (y_min+y_max)/2.0d0 )

	CALL moveto_w ( r0(1), ei0(1,j,k,l), wxy )
	do i=1+1,94
	result = lineto_w ( r0(i), ei0(i,j,k,l) )
c	status = setpixel_w( r0(i), ei0(i,j,k,l) )
	end do
	rgbcolor = #0000ff
	oldcolor = SETCOLORRGB( rgbcolor )
	CALL moveto_w ( r(1), ei(1,j,k,l), wxy )
	do i=1+1,118
	result = lineto_w ( r(i), ei(i,j,k,l) )
c	status = setpixel_w( r(i), ei(i,j,k,l) )
	end do
!
!     Print colored words on the screen.
!
	rgbcolor = #ffffff
	oldcolor = SETTEXTCOLORRGB( rgbcolor ) 
	CALL SETTEXTPOSITION( INT2(3), INT2(70), curpos )
	call outtext( ' The effect of i_coil .')
	write (6,fmt='(3x,3(a2,i2))') 'j=',j,'k=',k,'l=',l

	pause

	CALL CLEARSCREEN( $GCLEARSCREEN )
	end do
	end do
	end do

	CALL CLEARSCREEN( $GCLEARSCREEN )
	CALL SETVIEWPORT( INT2(30), INT2(30), xwidth-1-30, yheight-1-30 )
	CALL SETTEXTWINDOW( INT2(1), INT2(1), rows, cols )
	status = SETWINDOW( .true., x_min, y_min, x_max, y_max )
	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )

	do l=1,7
	do k=1,3
	do j=1,9
	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )
	status = RECTANGLE_W( $GBORDER, x_min, y_min, x_max, y_max )
	CALL MOVETO_W ( x_min, (y_min+y_max)/2.0d0, wxy )
	result = LINETO_W ( x_max, (y_min+y_max)/2.0d0 )

	call moveto_w ( r0(1), enl0(1,j,k,l), wxy )
	do i=1+1,94
	result = lineto_w ( r0(i), enl0(i,j,k,l) )
c	status = setpixel_w( r0(i), enl0(i,j,k,l) )
	end do
	rgbcolor = #0000ff
	oldcolor = SETCOLORRGB( rgbcolor )
	call moveto_w ( r(1), enl(1,j,k,l), wxy )
	do i=1+1,118
	result = lineto_w ( r(i), enl(i,j,k,l) )
c	status = setpixel_w( r(i), enl(i,j,k,l) )
	end do
!
!     Print colored words on the screen.
!
	rgbcolor = #ffffff
	oldcolor = SETTEXTCOLORRGB( rgbcolor ) 
	CALL SETTEXTPOSITION( INT2(3), INT2(70), curpos )
	call outtext( ' The effect of nl_coil .')
	write (6,fmt='(3x,3(a2,i2))') 'j=',j,'k=',k,'l=',l

	pause

	CALL CLEARSCREEN( $GCLEARSCREEN )
	end do
	end do
	end do

!
!     Print colored words on the screen.
!
	oldcolor = SETTEXTCOLORRGB( rgbcolor ) 
	CALL SETTEXTPOSITION( INT2(3), INT2(70), curpos )
	call outtext( ' The effect of m_coil .')

c	CALL SETVIEWPORT( INT2(0), INT2(0), xwidth-1, yheight-1 )
c	result=saveimage
c	&('fha.bmp',INT4(0),INT4(0),INT4(xwidth-1),INT4(yheight-1))

	return
	end

	subroutine paint_bf()
	use msflib
	implicit doubleprecision (a-h,o-z)
	dimension r(130),bf(120,4),bc(120,4)
	character fbf(4)*6
	common /cr/r/cbf/bf/cbc/bc

	LOGICAL               statusmode
	TYPE (windowconfig)   myscreen

	INTEGER(2)            status
	INTEGER(2)            xwidth, yheight, cols, rows
      
	INTEGER(4)            rgbcolor, oldcolor, rgbbkcolor, oldbkcolor
	TYPE (wxycoord)       wxy
	TYPE (rccoord)        curpos

	integer(4)            result

!     Set the screen to the best resolution and maximum number of 
!     available colors.
	myscreen.numxpixels   = -1
	myscreen.numypixels   = -1
	myscreen.numtextcols  = -1
	myscreen.numtextrows  = -1
	myscreen.numcolors    = -1
	myscreen.fontsize = -1
	myscreen.title = " "C
	statusmode = SETWINDOWCONFIG(myscreen)
	IF(.NOT. statusmode) statusmode = SETWINDOWCONFIG(myscreen)

	statusmode = GETWINDOWCONFIG( myscreen )
	xwidth  = myscreen.numxpixels
	yheight = myscreen.numypixels
	cols    = myscreen.numtextcols
	rows    = myscreen.numtextrows

c	rgbbkcolor = #ffffff
	rgbbkcolor = #000000
	oldbkcolor = SETBKCOLORRGB( rgbbkcolor )
	  
	CALL CLEARSCREEN( $GCLEARSCREEN )

c	CALL SETVIEWPORT( INT2(0), INT2(0), xwidth-1, yheight-1 )
	CALL SETVIEWPORT( INT2(30), INT2(30), xwidth-1-30, yheight-1-30 )
	CALL SETTEXTWINDOW( INT2(1), INT2(1), rows, cols )
	x_min=r(1)
	x_max=r(118)
	b_min=bf(1,1)
	b_max=bf(1,1)
	do n=1,4
	do i=1,118
	if (b_min.gt.bf(i,n)) b_min=bf(i,n)
	if (b_max.lt.bf(i,n)) b_max=bf(i,n)
	if (b_min.gt.bc(i,n)) b_min=bc(i,n)
	if (b_max.lt.bc(i,n)) b_max=bc(i,n)
	end do
	end do
	y_min=b_min-1.0d-1*(b_max-b_min)
	y_max=b_max+1.0d-1*(b_max-b_min)
	status = SETWINDOW( .true., x_min, y_min, x_max, y_max )
!     The 2.0_8 notation makes these constants REAL(8) 

c	rgbcolor = #000000
	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )  
	status = RECTANGLE_W( $GBORDER, x_min, y_min, x_max, y_max )
c     status = RECTANGLE( $GBORDER,INT2(0),INT2(0),xwidth-1,yheight-1 )
	fbf(1)='fbf_s1'
	fbf(2)='fbf_s2'
	fbf(3)='fbf_s3'
	fbf(4)='fbf_s4'
!
!     Plot the points.
!
	
	do n=1,4
	CALL SETVIEWPORT( INT2(30), INT2(30), xwidth-1-30, yheight-1-30 )
	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )  
	status = RECTANGLE_W( $GBORDER, x_min, y_min, x_max, y_max )
c	CALL MOVETO_W ( x_min, (y_min+y_max)/2.0d0, wxy )
c	result = LINETO_W ( x_max, (y_min+y_max)/2.0d0 )

	rgbcolor = #ffffff
	oldcolor = SETCOLORRGB( rgbcolor )  
	call moveto_w ( r(1), bf(1,n), wxy )
	do i=1+1,118
	result = lineto_w ( r(i), bf(i,n) )
c	status = setpixel_w( r(i), bf(i,n) )
	end do
	rgbcolor = #0000ff
	oldcolor = SETCOLORRGB( rgbcolor )  
	call moveto_w ( r(1), bc(1,n), wxy )
	do i=1+1,118
	result = lineto_w ( r(i), bc(i,n) )
c	status = setpixel_w( r(i), bc(i,n) )
	end do
!
!     Print colored words on the screen.
!
	rgbcolor = #ffffff
	oldcolor = SETTEXTCOLORRGB( rgbcolor ) 
	CALL SETTEXTPOSITION( INT2(5), INT2(60), curpos )
	call outtext( ' The bf. and bc. distribution .')
	write (6,fmt='(3x,a9,i2)') 'n_sector=',n

	CALL SETVIEWPORT( INT2(0), INT2(0), xwidth-1, yheight-1 )
	result=saveimage
	&(fbf(n)//'.bmp',INT4(0),INT4(0),INT4(xwidth-1),INT4(yheight-1))

	pause

	CALL CLEARSCREEN( $GCLEARSCREEN )
	end do

c	CALL SETVIEWPORT( INT2(0), INT2(0), xwidth-1, yheight-1 )
c	result=saveimage
c	&('fbf.bmp',INT4(0),INT4(0),INT4(xwidth-1),INT4(yheight-1))

	return
	end
